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Statement  of  the  Problem  Studied 

Increasing  the  service  lifetime  of  military  machinery  will  improve  logistics  in  the  field 
-(less  down-time  due  to  maintenance  or  replacement)  and  will  free  up  resources  for 
other  pressing  needs.  To  this  end,  the  Army  has  a  vested  interest  in  enhancing  the 
lifetime  of  gun  barrels,  which  are  limited  currently  by  erosion  of  the  steel  from  which 
the  gun  barrels  are  made.  However,  it  is  exceedingly  difficult  to  determine 
experimentally  in  situ  the  mechanisms  of  erosion.  Therefore,  computer  simulations 
can  help  elucidate  those  mechanisms  and  in  ultimately  extending  the  lifetime  of  gun 
barrels  by  discovering  new  ways  to  inhibit  erosion. 

As  a  result,  in  the  past  3  years  we  undertook  basic  research  to  create  new  simulation 
tools  and  to  investigate  interfacial  materials  chemistry  aspects  of  the  gun  tube 
erosion  problem  using  first  principles  atomistic  modeling  techniques.  In  brief,  we 
have  developed  minimal-adjustable-parameter  first  principles  models  of  iron, 
accurate  molecule-surface  interaction  theories,  large-scale  simulations  of  metals  via 
linear  scaling  methods,  and  chemical-mechanical  coupling  via  multiscale  modeling 
to  describe,  e.g.,  hydrogen  embrittlement  in  Fe.  In  terms  of  applications,  we  focused 
on  the  role  of  carburization:  we  examined  bulk  and  surface  properties  of  cementite, 
Fe3C,  thought  to  be  a  major  culprit  in  gun  tube  erosion  (GTE),  as  well  as  C 
incorporation  and  diffusion  in  Fe.  We  also  studied  ceramics  as  potential  alternative 
coatings  or  possibly  as  liners  to  protect  the  gun  tube. 

Summary  of  the  Most  Important  Results 

We  now  describe  the  most  important  findings  during  tenure  of  this  grant  from  our 
applications  of  existing  theoretical  techniques  to  problems  related  to  gun  tube 
erosion,  as  well  as  to  metal  particles  and  high  temperature  combustion  dynamics. 
Then  we  discuss  new  methods  we  developed  that  will  allow  more  realistic 
simulations  to  be  carried  but  in  the  future. 

First  Principles  Simulations  of  Carburization 

Gun  tube  erosion  involves  the  chemical  degradation  of  steel,  which  is  primarily  iron 
(doped  with  1%  carbon).  One  possible  mechanism  for  gun  tube  erosion  is 
carburization,  in  which  excess  carbon  is  incorporated  into  steel  (its  source  is  likely 
from  propellant  byproducts,  CO  and  C02).  Proposed  mechanisms  for  the  source  of 
carbon  include:  reduction  of  CO  by  hydrogen  to  produce  water  and  carbon  atoms, 
methane  dehydrogenation,  and  CO  disproportionation.  It  is  thought  that  the  carbon 
diffuses  into  the  steel  and  forms  cementite,  Fe3C,  which  has  a  lower  melting 
temperature  than  Fe  by  about  300  K.  Thus  one  form  of  erosion  is  thought  to  be  the 
formation  of  the  lower  melting  point  carbide,  which  then  liquefies  upon  gun  firing 
(temperatures  around  1700  K),  leading  to  erosion. 

We  studied  bulk  cementite  and  its  low-index  surfaces  using  first  principles  density 
functional  theory  (DFT).1  This  carbide  is  the  most  stable  of  all  iron  carbides,  but  it  is 
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still  metastable  with  respect  to  a-Fe  and  graphite.  Cementite  has  quite  a  complex 
crystal  structure,  with  16  atoms  in  the  orthorhombic  unit  cell:  four  6-coordinate  carbon 
atoms,  and  two  different  types  of  14-coordinate  Fe  atoms,  eight  of  one  type  and  four  of 
the  other.  Spin-polarized  DFT,  using  ultrasoft  pseudopotentials  and  the  generalized 
gradient  approximation  (GGA  -  critical  here  to  obtain  the  correct  bcc  ground  state  of 
bulk  Fe)  for  exchange  and  correlation,  yields  bulk  lattice  vectors  within  0.03  A  of 
experiment;  importantly,  the  lattice  expands  upon  incorporation  of  C,  which  will 
introduce  strain  at  the  interface  between  cementite  and  steel.  The  bulk  cohesive 
energy  is  within  6%  of  experiment.  We  predict  some  charge  transfer  from  Fe  to  C  in 
cementite,  but  it  is  still  metallic  (states  near  the  Fermi  level  are  Fe  d-states,  as  in  bulk 
iron).  We  also  find  a  reduced  magnetic  moment  compared  to  bulk  Fe  (this 
approximation  to  DFT  reproduces  the  Fe  magnetic  moment  to  within  0.25  jxB  and  the 
bulk  modulus  to  within  10%  of  experiment).2  The  surface  energies  of  Fe  (ranging  from 
2.3  -  2.6  J/m2)  follow  the  expected  trends  based  on  close-packing:  for  bcc,  the  (110)  is 
the  most  stable,  followed  by  (100),  and  then  (111)  is  the  least  stable;  for  fee,  the  most 
stable  is  the  (111),  followed  by  (110),  and  then  (100)  is  the  least  stable.  Seven  low- 
Miller-iridex,  stoichiometric  (non-polar)  surfaces  of  cementite  were  examined.  The 
relative  stability  of  these  cementite  surfaces  correlates  with  surface  smoothness:  the 
more  corrugated  the  surface,  the  higher  the  surface  energy  (the  latter  ranging 
between  2. 0-2. 5  J/m2).  The  surface  energies  of  both  Fe  and  cementite  are  very 
similar,  because  these  cementite  surfaces  are  dominated  by  Fe  termination.  All  but 
two  of  the  cementite  surfaces  have  surface  energies  lower  than  the  most  stable 
bcc(110)  surface  of  pure  iron.  This  suggests  that  lower  surface  energies  may  also 
provide  a  driving  force  for  formation  of  cementite  near  the  surface  of  steel,  leading  to  a 
surface  on  a  gun  tube  that  can  melt  at  lower  temperatures  than  steel  and  erode  away ,3 

Recently,  we  began  investigations  into  the  mechanism  of  carbon  incorporation  into 
iron.4  We  began  by  studying  C  incorporation  in  austenite  (fee  Fe),  since  the  solubility 
of  carbon  is  higher  in  austenite  than  in  ferrite  (bcc  Fe)  and  therefore  we  can  study  C 
incorporation  and  diffusion  in  austenite  with  smaller  simulation  cells.  We  first  verified 
that,  within  our  DFT  model  (all-electron  PAW  DFT-GGA  calculations),  it  is  exothermic 
to  dissolve  C  in  Fe.  We  found  that  it  is  exothermic  by  0.1  -  0.2  eV/atom,  depending  on 
the  concentration  (1-3%  C).  We  then  examined  the  diffusion  of  C  from  its  preferred 
octahedral  hollow  site  in  austenite  to  a  neighboring  octahedral  hollow.  We  find  a 
barrier  that  is  0.2  eV  too  low  compared  to  experiment  (1.26  eV  compared  to  1.47  eV). 
At  the  moment  that  is  at  a  -3%  concentration  level  (1  C  and  32  Fe  atoms  per  unit 
cell);  we  are  now  examining  the  barrier  at  a  lower  concentration  of  ~1%  (1  C  and  108 
Fe  atoms  per  unit  cell)  to  see  if  the  barrier  is  still  too  low.  This  provides  a  calibration  of 
the  DFT  exchange-correlation  energy,  so  that  we  can  know  approximately  what  errors 
we  will  expect  in  diffusion  barriers  when  we  begin  to  make  predictions  about  impurity 
diffusion  when  the  measured  values  are  not  available.  Next  we  will  study  C  diffusion 
in  bcc  Fe,  and  then  move  on  to  S  diffusion  predictions. 
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First  Principles  Exploration  of  Alternative  Protective  Coatings  for  Steel 

In  addition  to  examining  failure  mechanisms  for  gun  tube  erosion,  we  investigated 
possible  alternative  protective  coatings  for  gun  tubes.  Chrome  has  been  the  liner  of 
choice,  but  for  environmental,  as  well  as  other,  reasons,  alternative  protective 
coatings  are  desirable.  We  considered  ceramic  materials  as  candidates,  since 
ceramics  have  been  used  in  other  extreme  environment  conditions  quite  successfully 
to  protect  metal  components.  Ceramics  generally  exhibit  low  thermal  conductivity  and 
have  high  melting  points,  two  very  positive  traits  for  this  application.  In  addition,  we 
searched  for  ceramics  that  are  very  stable,  do  not  undergo  discontinuous  phase 
transitions,  have  similar  coefficients  of  thermal  expansion  to  that  of  Fe,  are  narrow- 
range  compounds  (no  variable  composition),  exhibit  good  mechanical  properties,  are 
environmentally  friendly,  and  potentially  adhere  strongly  to  iron.  After  scanning  known 
properties  of  various  ceramics,  we  selected  three  ceramics  on  which  to  focus:  TiC, 
ZrC,  (both  crystallize  in  the  NaCI  structure)  and  cubic  BN,  as  they  all  exhibit  excellent 
mechanical  properties  (hardness,  etc.),  have  very  high  melting  points  (>3000  °C), 
thermal  conductivities  at  least  one  order  of  magnitude  lower  than  that  of  the  current 
chrome  coating,  and  quite  similar  coefficients  of  thermal  expansion  to  that  of  steel  (~ 
10'5/°C),  DFT-GGA  theory  yields  bulk  properties  of  these  ceramics  (lattice  parameters, 
bulk  moduli,  and  energies  of  formation)  within  10%  of  experiment.5,6'7,8  Among  these 
crystals,  we  find  that  only  c-BN  exhibits  a  band  gap,  while  the  others  are  metals. 
Analysis  of  the  densities  of  states  indicates  that  TiC  and  ZrC  exhibit  a  mixture  of  ionic, 
covalent,  and  metallic  character,  while  c-BN  is  polar  covalent  in  nature.  We  observe 
the  effects  of  charge  transfer  in  our  studies  of  the  surfaces  of  these  ceramics,  where 
we  find  that  the  stoichiometric,  nonpolar  (100)  surfaces  of  TiC  and  ZrC  and  c-BN(IIO) 
are  by  far  the  most  stable.5,6,7,8  We  then  examined  the  adhesion  of  c-BN(IIO), 
TiC(100),  and  ZrC(100)  coatings  to  the  most  stable  face  of  a-Fe,  bcc  Fe(110).  The  bcc 
Cr-bcc  Fe  interface  was  also  investigated  for  comparison.  The  interfaces  for  the 
metal  carbide  coatings  with  Fe  are  dominated  by  strong  Fe-C  interactions,  while  the 
BN  coating  exhibits  stronger  Fe-N  than  Fe-B  interactions.  We  predict  that  TiC  is  the 
most  strongly  adhered  of  the  three  coatings,  with  an  ideal  work  of  adhesion  of  2.6 
J/m  .  ZrC  and  BN  coatings  have  much  weaker  adhesion,  -1.5  J/m2.  This  suggests 
TiC  is  the  best  candidate  for  an  alternative  to  Cr.  However,  even  the  TiC  work  of 
adhesion  is  lower  by  nearly  a  factor  of  2  compared  to  the  minimum  energy  required  to 
crack  Fe  (using  the  Griffith  criterion,  which  ignores  plastic  deformation  -  hence  the 
“minimum"  denotation),  which  suggests  that  the  coating  may  not  be  as  strongly 
adhered  as  might  be  needed.  Our  calculations  do  indicate  one  reason  why  Cr  is  a 
good  coating:  the  interfacial  adhesion  energy  of  Cr  to  Fe  (5.0  J/m2)  is  larger  than  the 
minimum  energy  required  (4.5  J/m2)  to  crack  Fe  (though  smaller  than  the  minimum 
energy  required  to  crack  Cr,  5.4  J/m2).  Though  the  iron-ceramic  bonding  is  weaker 
than  the  Cr-Fe  adhesion,  the  superior  protective  properties  of  these  ceramics  might 
warrant  them  to  be  considered  as  an  alternative  to  Cr. 
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First  Principles  Studies  of  the  Metal-lnsuiator  Transition  in  Metal  Nanoparticle 

Arrays 

Optica!  reflectivity  experiments9  showed  that  metal  nanoparticles  in  an  ordered  two- 
dimensional  array  undergo  a  transition  from  an  insulating  phase  to  a  conducting 
phase  upon  compression  of  the  array  (via  a  Langmuir  trough).  It  was  thought  that  the 
overlap  of  the  nanoparticle  wavefunctions  when  the  array  was  compressed  (leading  to 
smaller  particle-particle  separations)  led 'to  the  observed  electronic  transition.  Prior  to 
our  work,  first  principles  calculations  had  not  been  performed  to  try  to  determine  the 
mechanism  of  the  transition  because  the  number  of  atoms  required  to  simulate  such 
a  system  runs  into  the  many  hundreds  to  thousands  of  atoms.  We  have  a  DFT 
technique  that  can  treat  main  group  metals  within  a  linear  scaling  formalism,  as 
described  in  more  detail  later.  This  means  that  if  we  double  the  number  of  electrons 
treated,  we  only  double  the  cost  of  the  calculation.  Conventional  DFT  scales  cubically 
with  the  number  of  electrons  and  hence  only  small  (<100  atom)  particles  can  be 
treated. 

We  carried  out  a  series  of  calculations  on  a  hexagonally  periodic  array  of  AI249 
clusters,  with  two  such  clusters  in  the  unit  cell.  This  provides  nanoparticles  20  A  in 
diameter,  similar  to  the  size  of  the  particles  created  experimentally.  What  we  were 
able  to  do,  within  this  Hohenberg-Kohn  type  DFT,  is  to  self-conSistently  solve  directly 
for  the  density  of  the  -1500  valence  electrons  in  this  nanoparticle  unit  cell,  at  varying 
interparticle  separations.10  We  then  used  this  density  to  estimate  the  wavefunction 
overlap  between  each  nanoparticle,  and  this  was  in  turn  used  in  a  tight-binding  model 
to  ascertain  at  what  interparticle  separation  the  band  gap  closed.  We  find  a  large 
band  gap  at  large  separation  -  consistent  with  the  observed  insulating  behavior.  Then, 
as  the  separation  decreases  to  about  20%  greater  than  the  particle’s  diameter,  the 
band  gap  goes  to  zero  and  the  insulator-to-metal  transition  occurs.  This  is  entirely 
consistent  with  measurements  on  Ag  particles  of  similar  size  (26  A  diameter).  So  we 
have  captured  the  correct  physics  with  our  DFT  method.  Moreover,  we  noticed  that  the 
overlap  of  the  wavefunctions  at  the  point  at  which  the  gap  closes  is  extremely  small, 
which  means  that  very  tiny  overlaps  of  the  electronic  wavefunctions  between 
nanoparticles  are  sufficient  to  allow  conductivity  of  electrons  across  a  nanoparticle 
array. 

Ab  Initio  Molecular  Dynamics  of  Vinylidene-Acetylene  Isomerization 

As  a  first  step  toward  studying  the  dynamics  of  hot  gases,  we  examined  a  process 
important  in  acetylene  torches:  the  isomerization  of  vinylidene  to  acetylene  and  vice- 
versa.  It  turns  out  that  if  one  ignores  this  process  in  analyzing  the  kinetics  of  acetylene 
combustion,  one  derives  completely  the  wrong  thermochemistry  for  acetylene.  Thus 
this  reaction  is  an  important  one  for  the  complete  picture  of  acetylene  combustion. 
We  performed  complete  active  space  self-consistent  field  (CASSCF)  ab  initio 
molecular  dynamics  (AIMD)  simulations  of  the  preparation  of  vinylidene,  and  its 
subsequent,  highly  exothermic  isomerization  to  acetylene,  via  electron  removal  from 
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vinylidene  anion  (D2e=C~->  D2C=C:  ->  DC^CD).11  Remarkably,  we  find  that  only  20% 
of  the  vinylidenes  generated  this  way  isomerized,  suggesting  average  lifetimes  >  1  ps 
for  vibrationally  excited  vinylidene-d2.  Since  the  anion  and  neutral  vinylidene  are 
structurally  similar,  and  yet  extremely  different  structurally  from  the  isomerization 
transition  state  (TS),  neutral  vinylidene  is  not  formed  near  the  TS  so  that  it  must  live 
until  it  is  in  a  state  prepared  to. cross  the  barrier.  The  origin  of  the  delay  is  explained 
via  both  orbital  rearrangement  and  intramolecular  vibrational  energy  redistribution 
(IVR)  effects.  Unique  signatures  of  the  isomerization  dynamics  are  revealed  in  the 
anharmonic  vibrational  frequencies  extracted  from  the  AIMD,  which  should  be 
observable  by  ultrafast  vibrational  spectroscopy  and  in  fact  are  consistent  with 
currently  available  experimental  spectra.  Most  interestingly,  of  those  trajectories  that 
did  isomerize,  every  one  of  them  violated  conventional  transition  state  theory  by 
recrossing  back  to  vinylidene  multiple  times,  against  conventional  notions  that  expect 
highly  exothermic  reactions  to  be  irreversible.  The  dynamical  motion  responsible  for 
the  multiple  barrier  recrossings  involves  strong  mode-coupling  between  the 
vinylidene  CD2  rock  and  a  local  DCC  bend  mode  that  has  been  recently  observed 
experimentally.  Lastly,  the  multiple  barrier  recrossings  can  be  used,  via  a 
generalized  definition  of  lifetime,  to  reconcile  extremely  disparate  experimental 
estimates  of  vinyiidene’s  lifetime  (differing  by  at  least  six  orders  of  magnitude).  Thus, 
we  once  again  demonstrated  the  usefulness  of  ab  initio  molecular  dynamics  to  reveal 
behavior  far  beyond  what  a  simple  study  of  the  potential  energy  surface  would  predict. 
Such  studies  will  be  carried  out  for  the  hot  propellant  gases  involved  in  gun  tube 
erosion  in  the  future. 

Development  of  Electronic  Structure  Theories  for  the  Condensed  Phase 

This  section  describes  the  continued  development  of  good  solid-state  models  for 
metals  and  molecules  interacting  with  metals.  This  work  consists  of  improvements 
to  both  the  scaling  of  density  functional  theory  (DFT)12  via  the  development  of  accurate 
kinetic  energy  functionals,  3,14,15  the  accuracy  of  ground  state  DFT  via  systematic 
corrections  to  it  in  a  local  region  (e.g.,  where  a  chemical  bond  is  breaking  or  forming), 
as  well  as  the  ability  to  calculate  local  excited  states  on  or  in  metals. 16,17,18  It  also 
consists  of  improved  descriptions  of  the  effect  of  the  core  electrons  on  the  valence 
electrons  via  spin-dependent  pseudopotentials ,19'20  We  are  also  working  linking  our 
atomistic  simulations  to  larger  length  scale  mechanical  engineering  simulations  of 
cracking. 21,22,23  On-the-fly  linking  of  our  linear  scaling  DFT  theory  for  metals  to  the 
quasicontinuum  method  for  treating  deformations  of  solids  has  been  achieved.24 
This  will  allow  more  realistic  simulations  of  large-scale  gun  tube  erosion  to  be 
tackled  in  the  future. 

Kinetic  Energy  Density  Functionals  (KEDF’s) 

It  is  possible  to  approximately  solve  the  Schrodinger  equation  directly  for  the  electron 
density.  In  order  to  do  this,  one  must  have  not  only  an  expression  for  the  exchange- 
correlation  energy,  but  also  one  for  the  kinetic  energy,  in  terms  of  the  density. 
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Recently,  we  developed  new  KEDF’s  that  satisfy  both  smoothly  varying  and  rapidly 
varying  density  limits  of  exact  linear  response  theory.13,14  These  new  functionals  yield 
accurate  structures,  energy  differences,  vacancy  formation  energies' ‘densities,  and. 
even  surface  energies  and  densities  (where  accuracy  is  relative  to  Kohn-Sham'  DFT, 
since  the  KEDF  s  are  an  approximation  to  K-S  DFT).  The  advantage  of  the  density-only1 
DFT  method  is  that  the  cost  of  the  method  scales  linearly  with  the  size  of  the  system, 
allowing  one  to  solve  for  the  electron  densities  and  properties  of  thousands  of  atoms; 
moreover,  it  is  the  only  such  method  available  for  treatment  of  metallic  systems.  As 
mentioned  above,  we  used  these  new  functionals  to  examine  the  electronic  structure 
of  nanoparticle  arrays,  work  that  would  have  been  impossible  with  K-S  DFT  and  would 
not  have  been  accurate  with  previous  KEDF’s.13  Further  work  is  in  progress  to  extend 
these  KEDF’s  to  nonmetallic  systems,  but  this  is  an  extremely  difficult  problem. 

We  recently  wrote  a  review  of  all  the  issues  associated  with  KE  DFT,15  where  we  make 
it  clear  that  the  two  outstanding  issues  are:  (i)  how  to  constrain  these  KE  functionals 
further  to  treat  localized  electron  densities  and  (ii)  accurate  pseudopotentials  (which 
represent  the  screened  ion-valence  electron  interaction)  need  to  be  developed  that 
are  local  in  nature  (i.e.,  only  depend  on  one  coordinate  and  do  not  involve  projection 
operators  that  must  act  on  a  wavefunction).  Regarding  the  first  point,  we  are 
examining  the  use  of  the  idempotency  constraint  on  the  density  to  see  if  this  will  help 
constrain  the  KE  functionals,  to  make  them  applicable  to  nonmetallic  systems.  Our 
efforts  in  the  last  two  years  have  focused  more  intensely  on  the  second  point, 
because  the  way  to  proceed  was  clearer.  We  developed  local  pseudopotentials 
derived  by  inverting  the  Kohn-Sham  equations  and  finding  a  pseudopotential  that 
reproduces  a  Kohn-Sham  (orbital-based  DFT)  density  for  both  atoms25  and  bulk 
environments.26  In  particular,  we  have  shown  that  accurate,  transferable  local 
pseudopotentials  for  an  element  that  is  not  at  all  nearly-free-electron-like  can  be 
developed  by  inversion  of  the  Kohn-Sham  equations  and  exploiting  the  first 
Hohenberg-Kohn  theorem.  The  latter  states  that  given  a  density  there  is  a  one-to-one 
mapping  of  the  exact  potential  associated  with  that  density.  Thus,  given  a  high  quality 
density  (we  use  Kohn-Sham  DFT  densities  from  nonlocal  norm-conserving 
pseudopotential  calculations),  one  should  be  able  to  find  the  exact  Kohn-Sham 
effective  potential  that  reproduces  that  density.  Unfortunately,  this  is  an  inverse 
problem  with  no  unique  solution.  We  tried  many  algorithms  that  have  been 
suggested,  and  found  the  one  due  to  Wang  and  Parr  to  be  the  most  stable.  We  first 
utilized  this  approach,  which  updates  a  guessed  potential  utilizing  the  target  density, 
the  current  iteration's  orbitals  and  eigenvalues,  in  the  context  of  atoms  only.  We  then 
worked  out  a  scheme  to  do  the  same  thing  given  a  target  bulk  density  (this  is 
nontrivial  as  the  only  data  points  available  are  at  the  Bragg  vectors  for  a  given  crystal 
structure).  We  now  have  a  robust  algorithm  in  the  solid  state  for  this  purpose  and 
have  shown  that  we  can  produce  accurate,  transferable  local  pseudopotentials  that 
reproduce  nicely,  within  Kohn-Sham  DFT,  the  phase  orderings  of  8  different  phases  of 
Si  predicted  by  nonlocal  pseudopotentials  within  Kohn-Sham  DFT.  When  we  now 
use  this  local  pseudopotential  in  orbital-free  DFT  calculations  of  different  phases  of 
Si,  any  errors  incurred  have  now  been  isolated  to  only  the  kinetic  energy  density 
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functional.  This  provides  the  means  to  then  move  forward  jn  improving  KE  density 

functionals  for  nonmetals. 

*  -  '  : 

As  a  last  note,  in  the  summer  of  2002,  there  was  workshop  at  CECAM  in  Lyon,  France, 
devoted  to  kinetic  energy  density  functionals.  The  P.l.  was  invited  to  speak,  but  instead 
sent  her  former  postdoc,  Alex  Wang,  and  a  current  graduate  student,  Baojing  Zhou  to 
the  meeting.  Alex  Wang  talked  about  these  new  kinetic  energy  density  functionals  and 
their  applications,  and  Baojing  Zhou  presented  a  poster  on  deriving  ab  initio  local 
pseudopotentials  based  on  Kohn-Sham  nonlocal  pseudopotential  bulk  densities. 
The  P.l.  later  heard  from  more  than  one  senior  participant  that  the  P.I.’s  KEDF’s 
remain  the  state-of-the-art  and  the  most  advanced/accurate  available  today. 

Spin-Dependent  Pseudopotentials 

In  order  to  correctly  model  Fe,  Cr,  and  compounds  of  Fe,  it  is  critical  to  capture  the 
ferromagnetic  nature  of  Fe  and  (at  least)  the  local  anti-ferromagnetic  behavior  of  Cr 
correctly  within  DFT,  since  DFT  is  our  starting  point  for  all  other  calculations.  Efficient 
solid-state  DFT  calculations  employ  plane-wave  basis  sets,  which  describe  the 
valence  electrons,  and  pseudopotentials,  which  represent  the  nucleus  and 
associated  core  electrons.  Generally,  a  nonmagnetic  reference  state  is  used  to 
construct  such  pseudopotentials.  Unfortunately,  when  such  pseudopotentials  are 
used  for  transition  metals,  large  errors  occur  in  the  description  of  the  d-electrons 
(though  the  nonlinear  core  correction  to  exchange-correlation  reduces  these  errors,  at 
the  expense  of  introducing  more  parameters).  We  have  ameliorated  this  problem  by 
developing  spin-dependent,  norm-conserving  .  pseudopotentials.27  These 

pseudopotentials,  containing  both  a  nonmagnetic  term  and  a  part  that  depends  self- 
consistently  on  the  spatially  varying  net  spin  polarization  (with  no  additional 
parameters),  greatly  improve  the  transferability  of  the  pseudopotential  compared  to 
traditional  formulations.  Not  only  do  they  accurately  reproduce  atomic  state  splittings 
for  the  entire  first  and  second  row  transition  series,  but  we  have  shown  they 
reproduce  reasonably  well  the  structural  properties  of  bulk  Fe,  Co,  and  Ni.19  Moreover, 
the  predicted  magnetic  moment  is  much  closer  to  experiment  when-  the  spin- 
dependent  pseudopotentials  are  employed.  We  have  now  implemented  these 
pseudopotentials  in  a  reciprocal-space  bulk  DFT  code  (CASTER),28  which  is 
-  considerably  more  efficient  than  the  real  space  bulk  code  in  which  we  first 
implemented  the  method.  The  advantage  of  these  pseudopotentials  is  that  they 
provide  the  same  degree  of  accuracy  without  introducing  a  plethora  of  empirical 
parameters  needed  by  other  approaches.  The  disadvantage  is  they  have  a 
somewhat  higher  cost  associated  with  their  use.  To  the  extent  that  one  prefers  to 
minimize  the  number  of  empirical  parameters  in  a  method,  this  approach  is 
aesthetically  more  pleasing.  In  addition,  the  embedding  theory  below  requires  the 
use  of  so-called  norm-conserving  pseudopotentials,  which  necessitates  ,  the  use  of 
these  (norm-conserving)  spin-dependent  pseudopotentials  for  Fe-/;ontaining 
materials.  - 
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Embedding  Theory  for  Accurate  Energetics  in  a  Local  Region  of  Condensed 

Matter 

While  one  can  go  far  studying  metals  using  the  state-of-the-art  in  condensed  matter 
electronic  structure,  namely  periodic  DFT,  there  are  times  when  DFT  fails  inexplicably, 
due  to  the  approximations  made  in  accounting  for  electron  exchange  and  correlation! 
We  developed  a  method  that  improves  upon  the  DFT  result  by  means  of  explicit 
convergence  of  electron  correlation  (via  perturbation  theory  and  Cl).  The  Army 
Research  Office  (ARO)  helped  support  the  development  of  this  theory  where  a  metal 
cluster  is  embedded  into  a  periodic  slab  representation  of  the  infinite  metallic 
conductor.  We  demonstrated  that  this  systematically  improves  the  description  of 
electron  correlation  in  a  local  region,  over  that  offered  by  DFT.29  In  practice,  a  local 
region  of  the  solid  is  represented  by  a  small  cluster,  for  which  high  quality  quantum 
chemical  calculations  are  performed.  The  interaction  of  the  cluster  with  the  extended 
condensed  phase  is  taken  into  account  by  an  effective  embedding  potential.  This 
potential  is  constructed  within  periodic  DFT  and  is  used  as  a  one-electroh  operator  in 
subsequent  cluster  calculations.  This  was  done  originally  within  the  context  of  HF  and 
Maller-Plesset  perturbation  theories  (MP-n)  through  fourth  order  for  improved 
descriptions  of  electron  exchange  and  correlation  in  electronic  ground  states,  but 
recently  we  extended  it  to  provide  also  a  fairly  accurate  description  of  localized  excited 
states  in  condensed  matter. 16,17,18  Such  local  excited  states  have  not  been  accessible 
from  ab  initio  theory  for  metallic  systems  before  this  work.  Among  a  variety  of 
benchmark  calculations,  we  investigated  a  CO  molecule  adsorbed  on  a  Pd(1 1 1) 
surface.  By  performing  CASSCF,  Cl,  and  MP-n  theories,  we  not  only  were  able  to 
obtain  accurate  adsorption  energies  via  local  corrections  to  DFT,  but  also  vertical 
excitation  energies,  e.g.,  for  an  internal  (5o— >27r*)  excitation  within  the  adsorbed  CO 
molecule.  We  have  shown  that  the  agreement  with  experimental  electron  energy  loss 
spectra  (EELS)  improves  dramatically  upon  use  of  the  Cl  embedding  theory, 
compared  to  Cl  calculations  on  finite  clusters.  Thus,  we  demonstrated  that  our  new 
scheme  provides  fairly  accurate  predictions  of  local  excited  states  on  metal  surfaces, 
as  well  as  a  systematic  means  of  improving  locally  on  ground  state  properties. 

Work  in  progress  on  the  embedding  theory  involves:  (i)  calculation  of  forces  such  that 
CASSCF  dynamics  of  molecules  impinging  and  reacting  on  metal  surfaces  will  be 
possible  (the  coding  for  this  is  nearly  done  and  partially  tested);  (ii)  exploration  of  the 
robustness  of  the  kinetic  energy  functional  potentials  used  in  the  embedding 
potential,  when  applied  to  other  materials  (currently  we  are  examining  this  for  Al,  Cu, 
and  Pt  surfaces);  and  (iii)  extension  of  the  embedding  theory  to  MRSDCI  (currently  we 
have  developed  the  ability  to  get  Cl  natural  orbitals  and  are  implementing  the 
calculation  of  Cl  densities  from  the  natural  orbitals,  so  that  we  can  construct  the 
cluster  only  terms  of  the  embedding  potential  from  a  MRSDCI  density). 

This  theory  will  be  used  to  refine  diffusion  and  reaction  energetics  on  Fe  surfaces  and 
to  examine  propellant  by-product  surface  and  bulk  chemistry  on  Fe  in  the  years  to 
come. 
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Bridging  Length  Scales  in  Materials  Simulations 

•  We' are  funded  elsewhere  to  work  on  the  bridging  length  scales  problem,  however  the 
work  in  that  arena  will  have  considerable  impact  on  our  ability  to  model  gun  tube 
erosion,  and  so  our  progress'  in  this  area  is  briefly  described  here.  First,  we  are 
collaborating  with  Michael  Ortiz,  a  Professor  in  the  Graduate  School  of  Aeronautics  at 
Caltech,  to  link  our  atomistic  information  to  micron-to-mm  scale  materials  failure. 
This  is  being  tackled  in  two  ways.  The  first  is  by  what  may  be  referred  to  as  an 
informed  continuum  model,  in  which  information  from  the  atomic  scale  is  used  as 
input  into  a  separate  continuum  level  simulation  of  crack  propagation.  We  have 
calculated  the  force  laws  that  describe  how  cracks  form  in  metals  and  oxidized  or 
embrittled  metals  and  have  formulated,  via  renormalization  group  theory,  a  means  of 
rescaling  the  atomic  scale  forces  involved  in  crack  opening  up  to  an  effective  force  law 
for  the  microscale.  This  work  promises  to  introduce  chemical  effects  on  mechanical 
materials  failure  for  the  first  time.23  Work  in  progress  includes  calculating  the  critical 
opening  displacement,  maximum  tensile  stress,  and  work  of  adhesion  as  a  function 
of  hydrogen  concentration  in  steel.  We  are  developing  a  coupled  diffusion/cracking 
model  for  hydrogen  embrittlement  in  steel.30 

Secondly,  we  developed  an  on-the-fly  coupling  of  electrons/atoms  to  a  micron  scale 
mechanical  engineering  model,  via  linking  our  linear  scaling  orbital-free  DFT  code  to 
a  local  quasicontinuum  code  (which  seamlessly  goes  from  a  finite  element  averaged 
description  of  millions  of  atoms  down  to  nanometer  resolution  where  needed).  What 
is  required  is  the  ability  to  treat  bulk  unit  cells  of  arbitrary  distortions  and  to  compute 
the  stress  tensor  on  that  cell,  neither  of  which  was  implemented  in  our  OF-DFT  code. 
The  idea  is  that  the  deformations  of  each  DFT  cell  are  determined  by  the 
quasicontinuum  method  while  the  DFT  calculations  provide  the  stress  tensors  to  the 
quasicontinuum  code  to  determine  how  to  next  deform  each  finite  element  node.  We 
have  carried  out  nanoindentation  simulations  of  A!  thus  far.24  We  should  be  able  to 
predict  hardness  of  metal  alloys,  in  the  presence  of,  e.g.  precipitates,  in  this  way. 

Thirdly,  we  are  developing  an  informed  continuum  approach  to  martensitic  phase 
transformations  in  steel,  by  calculating  the  elastic  properties,  interfacial  energies,  and 
transition  paths  for  pressure-induced  bcc  to  hep  transitions  in  Fe.31 

These  new  software  tools  will  allow  us  to  model,  e.g.,  shock  wave  deformations  in  Fe 
and  chemically-induced  cracking  of  Fe. 
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Inventions 

Development  of  the  following  algorithms  and  codes: 

•  Local  pseudopotentials  for  use  in  our  orbital-free  density  functional  method, 
derived  by  inverting  the  Kohn-Sham  equation 

•  Generalization  of  our  orbital-free  density  functional  code  to  treat  non-orthorhombic 
cells  and  to  compute  the  stress  tensor,  in  order  that  arbitrarily  deformed 
simulation  cells  may  be  treated 

•  spin-dependent  pseudopotentials  in  the  Kleinman-Bylander  form  for  use  in  both 
real  and  reciprocal  space  solid  state  density  functional  theory  codes. 

•  An  embedding  theory  of  ground  and  excited  states  in  bulk  metals  and  metal 
surfaces 

•  Coupled  orbital-free  density  functional/finite  element  code  for  multiscafe  modeling 
of  deformations  of  micron  size  samples  of  metal 
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